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A novel strategy is introduced in order to include variations of the nonlinearity into the nonlinear 
schrodinger Equation. This technique, which relies on renormalization, is in particular well adapted 
to nanostructured optical systems where the nonlinearity exhibits large variations up to two orders 
of magnitude larger than in bulk material. We show that it takes into account in a simple and 
efficient way the specificity of the nonlinearity in nanostructures that is determined by geometrical 
parameters like the effective mode area and the group index. The renormalization of the nonlinear 
schrodinger Equation is the occasion for physics oriented considerations and unveils the potential of 
Photonic Crystal waveguides for the study of new nonlinear propagation phenomena. 
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I. INTRODUCTION 

Nanostructured optical systems offer a great plastic¬ 
ity because most of their optical properties can be engi¬ 
neered through the design 00 . Thus the control by 
the geometry allows one to get on demand properties, 
precisely adapted to the phenomena one wants to inves¬ 
tigate p, I Til 1 l(il . This is both interesting for applications 
where it is then possible to optimize all optical functions 
to a large extend EH and for fundamental investigations 
where the experimental conditions could be precisely set 
to maximize a given effect. 

This plasticity comes usually along with large variations 
of the optical properties over a small bandwidth because 
effects of geometry depend highly on the actual wave¬ 
length size. For instance in Photonic Crystal waveg¬ 
uides, variations of the nonlinear effective Kerr effect up 
to 100% can occur over a lOnm bandwidth [T^] . Con¬ 
sequently the modeling of the nonlinear propagation of 
optical pulses in such structures is essential for under¬ 
standing precisely the interplay between the different ef¬ 
fects taking place; such a task is challenging 0 - 10 ]. 

The use of the generalized Nonlinear Schrodinger Equa¬ 
tion (GNLSE) is usually preferred to direct nonlinear 
FDTD simulation, though that latter being more accu¬ 
rate [H], [23| , because it has very low computation bur¬ 
den and also provides a direct link between the effective 
coefficient in the GNLSE and the phenomena that are 
observed. Namely it is straightforward to add phenom¬ 
ena specifically found in semiconductor optics, like for in¬ 
stance the effect of nonlinear absorption and free carriers 
[Til 0 ] . Moreover discussion and interpretation are made 
easier as the coefficients in the GNLSE are derived into 
effective lengths like the dispersion length Ld = T§ / fc, 
the nonlinear length Ljssl = 1/(7 Po) or the shock forma¬ 
tion length L sh ock = 0.43Latl/|tatl| 0, etc. The effec¬ 
tive lengths [ 29 } provide a rapid overview on the relative 
strength of the different competing effects. However, be¬ 
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cause the GNLSE relies on several approximations, it is 
important to check that the physics governing the optical 
system is still correctly described. 

Higher order nonlinear effects, for weak perturbations, 
are taken into account by adding successive corrective 
terms to the initial NLSE. Here we focus on the inclusion 
of variations of the nonlinear Kerr response with the an¬ 
gular frequency. Using the transformation Aw —> *<9j[30|], 
the Taylor expansion of the Kerr coefficient (hence its 
variation in the frequency domain) can be included into 
the GNLSE (time domain). Such decomposition -limited 
to the first order- takes into account self-steepening ef¬ 
fects and can be used to model the formation of optical 
shock front 0, |3l|. 

Despite being easy to implement, such technique suf¬ 
fers from two major limitations. 

First it uses the derivation of an aggregated effective 
parameter computed for a given frequency, whereas by 
essence nonlinearity involves also the interactions be¬ 
tween several waves at different frequencies. Notably, not 
only the effective mode area, but also the variations of 
the modes overlap -directly related to cross phase mod¬ 
ulation (XPM)- are important in integrated optics. A 
simple Taylor expansion of the aggregated Kerr coeffi¬ 
cient cannot take into account that latter effect. This is¬ 
sue was addressed in Ref. 0] where it was demonstrated 
that the physics related to variations of the mode area 
overlap could be in principle retrieved by the addition of 
extra operators; and that such corrections could actually 
have a noticeable impact on pulse propagation. 

The second limitation is due to the nature of the Taylor 
expansion itself: it is intended for weakly varying func¬ 
tions and convergence tends to be slow as soon as the 
functions exhibit large or non-monotonic variations. In 
case of large variations of the nonlinear effective coef¬ 
ficients, the best option would then be to split the to¬ 
tal bandwidth into sub-domains wherein variations of 
the nonlinear properties are negligible 0, 0. Thus 
each sub-domain is described with good accuracy us¬ 
ing a GNLSE; then the different equations are connected 
through ad-hoc parameters that describe with precision 
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effects like cross phase modulation (XPM) or four wave 
mixing (FWM). This method is especially well suited 
to describe FWM because for most practical applica¬ 
tion pumps, signal and idler spectra are well separated 
; moreover the small number of beams involved allow 
setting up a set of coupled mode equations with only a 
few equations. However, if the spectrum is continuous 
over a large bandwidth, splitting the initial simulation 
domain into several pieces is complex because there is no 
obvious choice for the separation frequency between the 
sub-domains; and numerical artefact might appear at the 
junction between two sub-domains. 

Consequently neither of the two methods gives sat¬ 
isfactory results in case the single pulse bandwidth 
extends over large variations of the nonlinear effec¬ 
tive parameters. This is especially true for dispersion- 
engineered Photonic Crystal waveguides (PhCWGs) [T|| 
where the nonlinearity arises mainly from the slow light 
enhancement (l9j] . Indeed any changes of the group in¬ 
dex is immediately reported through the slow light fac¬ 
tor S 2 = (rip / tt-o ) 2 [341 on the Kerr nonlinearity, which 
is expressed as 7 e // = u)ri 2 i/(cA e ff) * S 2 . Variations 
of the modal area are also important in PhCWGs, but 
still much weaker than those associated with the slow 
light variation. Note that while S ranges from 2 to 10 
in PhCWGs [24[, its effect is less sensitive in Nanowires 
on silicon which in contrast are much more impacted by 
variations of the modal area. Indeed the slow light effect 
in Nanowires is roughly S' ~ 1.2 [35} and is in fact al¬ 
ready implicitly taken into account in the expression of 
the effective area jTTot. 

The method we present here consists in defining and 
solving the GNLSE for another field than the power flux 
(cf. Fig |T]). The point is to find the pseudo held that 
leads to the GNLSE exhibiting minimal variation of the 
nonlinear coefficient. Especially this technique is well- 
suited for nanostructured optical systems where the non¬ 
linearity is defined by the effective area and the group in¬ 
dex. Previous studies focused on the ab initio definition 
of the effective parameters of the propagating equation 
once the eigen modes profile is known. It was demon¬ 
strated that the presence of slow light must be taken 
into account in the normalization in order to guarantee 
that the correct energy how through the medium is pre¬ 
served despite the presence of strong material dispersion 
or slow light [37j. This hrst normalization step is cru¬ 
cial. However, the resulting nonlinear propagation equa¬ 
tion will contain also a few pre-factors; this aspect is not 
important for few modes problems like FWM or para¬ 
metric generation where one can dehne a set of coupled 
equations, but becomes problematic for continuous mode 
problems. For example the propagation equation can be 
turned back into the form of a classical NLSE, but at the 
cost of extra approximations which ultimately alter the 
accuracy of the numerical solution. The second step con¬ 
sists then in the renormalization of the energy how into a 
variable whose nonlinear propagating equation has a sim¬ 


pler form. Our hnding is that instead of computing the 
propagation of the power flux, it is better to weight it by 
the photonic density of states. The group index plays an 
important role in this second normalization step because 
the hrst derivative of the band diagram corresponds at 
the same time to the group velocity and to the optical 
density of states. 

First, we will hrst review how the slow light effect is intro¬ 
duced into the nonlinear propagation equations; and thus 
insist on the two contributions of slow light on the non¬ 
linear enhancement. We will then introduce our renor¬ 
malization technique and apply it to the computation 
of nonlinear pulse propagation in dispersion-engineered 
PhCWGs. 


II. NLSE 

We start from the Maxwell’s equations. As we are deal¬ 
ing with unidirectional guided propagation (i.e. waveg¬ 
uides), we have split the spatial coordinates into trans¬ 
verse coordinates f± and z along the waveguide direction. 

V x E(r±,z,t) = -Mo d t H(fl,z,t) (1) 

V x H(r±,z, t) = e(r±,z ) d t E(r±,z,t ) + <9 t P/v L (E) 

( 2 ) 

When Pnl = 0, the above equations accept a set of 
eigen solutions in the frequency domain: 

e(rl, 2 ,w) = e b ioch(rj_, z - n.a,uj)e~ luJt+lK( ' ul)z (3) 

h(r±,z, u) = h b i och (r _l, z - n.a , u )e~ wb * K M z (4) 

with n G Z and a being the PhC lattice parameters. 
K(to) corresponds to the propagation constant at the fre¬ 
quency u>. For periodic waveguides like PhCWGs, the 
electric and magnetic field distribution is determined by 
a Bloch mode which is a periodic function with periodic¬ 
ity a. Hereafter, the dependence on (r±,z) of the eigen 
and Bloch modes is implicitly assumed and we will use 
the short notation e(u>), h(u>), e B ioch(v), h B i och (oj). 
Similarly as detailed in [38|, the forward propagating 
equation is obtained by multiplying eqs. m by the solu¬ 
tion of the unperturbed ( Pnl = 0) system (i.e. e(w)* and 
h(u>)*). Substracting to each other the two new equations 
leads, after few algebra -and especially using the identity 
AV x B = V.(A x B) + B.V x A-, to: 

V.(E x h(u)* — H x e(w)*) = 

-e(u>)d t pNL{E) - d t (e 0 eE.e(oj) + H.h(u)) (5) 

We will now look for a perturbative solution, in the 
frequency domain, of the form 
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E(rl, z, Aw) = A(z, Au)e-' A ^S B i oeh {u)e-' u ”t t+ ' K t u "f* ( 6 ) 

z, Aw) = A(z, (7) 


Where AAT(w) = K{w) — K(io re f ) and Aw = w — 
uj re f. Note that the phazor term in eqs. m is different 
from the one found in eqs. m , hence we would like 
to compute directly the evolution of the pulse envelope 


centered at frequency u re /. Also we assume the envelope 
d z A/A « 1/a varies slowly, so d z A(z,u>) is constant 
over a single PhC lattice. We inject eqs. IbllTl) into eq. [5] 
and integrate over a unit PhC cell; this leads to: 


(d z A(z, Aw) — iAK(uj)A(z, Aw)) 


^ Bloch (cj) X h]3ioch(bj) h/Blochi^) ^ ^Bloch{^))dv 


cell 


•JJL 


eBioch(u)* ■ Pnl{Au>, E)dr 3 


( 8 ) 


A(z, Aw) is defined with regard to a central frequency 
reference w re /. Thus eq. [S] gives directly access to the en¬ 
velope of the pulse, without the fast oscillations in space 
and time {w re /, K(ui re f)}. The integral in the left hand 
side of eq. [ 8 ] corresponds to twice the integral of the 
Poynting vector II W over the PhC cell. The link between 
the Poynting vector and the Bloch mode energy is set 
through the relationship: 


W w is a normalization factor, which corresponds 
to the Bloch mode energy l/2fJJ ceU e\e B i och (uj)\ 2 + 
fi\h B i och (uj)\ 2 dr 3 . Finally we have 



x h^dr 3 = 


z ■ I Pdr 2 * a 


Surface 
Vg(U})W U 


(9) 


J 


d z A(z, Aw) 


,akma(z,a u ) + , 1 ^L w ^ 



e B ioch{u) * .PNL(Au>,E)dr 3 


( 10 ) 


Before going further into details, let us focus an instant 
on the peculiar normalization choice for both A(z, w) and 
W u . In the linear regime, the choice is to preserve the 
power flux, so the computed field directly corresponds to 
the energy flowing though the medium. Using eq. [9] one 
finds that the total energy transiting through the waveg¬ 
uide at frequency w is P(z,u) = \A(z, w)| 2 v s (w)W w /a. 
As a result the Bloch modes are normalized such as 
Vg{uj)W u /a = 1; so one gets directly from the NLSE 
P = \A\ 2 . However such normalization is done indepen¬ 
dently of the nonlinear problem which is considered; con¬ 
sequently it may not necessary be the best choice from a 
pure numerical point of view. 


Regarding the nonlinear polarization Pnl, we as¬ 
sume that it is due to the x response ex¬ 
pressed as: P NL {r , wo) = 3/2eoXmj wl _ 2 _ 3 (P*(r-)u, 1 • 
E(r) ul2 ).E(r) ul3 )S(u }2 + W 3 — wi — wo) , where <5 stands 
for the Dirac delta function. The exact form of Pnl is 
directly related to the nonlinear tensor Q; and conse¬ 
quently will have a different formulation depending on 
the material that is actually considered. However the 
conclusions presented here are general and could be eas¬ 
ily extended to any peculiar form of Pnl- Finally, we 
obtain 
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d z A(z,u > 0 ) = iiC(wo)-4(wo) 


lU}ori2l \/ n guo 
cA e ff{uj 0 ) nl 


E(z,uJi)*E(z,u}2)E(z,u] 3 )gh ul0 (uJi,u]2,U3)8( u} 2 + w 3 - uq - ui 0 )duif 


( 11 ) 


We have defined n^i = 3\l 3 V(4ceon§). We also intro¬ 
duced the error function gh u which takes into account 
the fact that the mode overlap integrals may deviate from 
A e ff(oJo) for frequency uq different than uiq. Hence 


A eff(u o) 


{IffCell e ( r )l e ^o| 2 ^ r3 ) 2 

an tfffNLmat .\ e ^\ Adr3 


( 12 ) 


gh u ,o(u!i,U!2,U!3) = 


fffNLmatf^O ' ' e u 3 )dr 3 


Iff] 


NLmat, 


|e^ 0 | 4 dr 3 


(13) 


NLzmat indicates that the integral is performed over 
the nonlinear material in the PhC unit cell. At this point, 
one notes that the nonlinearity depends on the intensity 
of electric field E while in eq. [TT] the nonlinear pulse 
evolution is set for the power flux P = \A\ 2 |u|. Eq. [9] 
is then used once more and one finally gets: 


d~A(z,w 0 ) = iK(uj 0 )A(uj 0 ) + i 


u) 0 n 2 i y/n^L 

cAeffiu 0 ) nl 


A(z, uii)*A(z, uj 2 )A(z , u 3 )gh u0 {u i, w 2 , w 3 ) 

Jn g {ui)n g (u>2)n g (u} 3 )8{u2 + u) 3 -u >i - uJo)du 3 (14) 


Unfortunately, the ^n g (u)i= 1 , 2 , 3 ) and ghu,o(uJi, ^ 2 , w 3 ) 
terms in the right hand side depend on <^= 1 , 2 , 3 , hence 
these prefactors cannot be set aside the integral. How¬ 


ever, the formulation of a standard NLSE equation - 
expressed in the frequency domain- can be retrieved as¬ 
suming that gh u , = 1 and that ^Jrlg(uJ^) do no vary: 


d z A(z,uj 0 ) = iK(lo 0 )A(z,u} 0 ) + 1 a; ° rt2J f A(z,u} 1 )* A(z,U} 2 )A(z, v 3 )S(uj 2 +u 3 - uq - u} 0 )du 3 (15) 

cAeff(uj 0 ) nl J 


Note that eq. [15] could have also been directly ob¬ 
tained through the derivation of eq. [TT] for a monochro¬ 
matic wave propagation. In fact this is most of the time 
how the effective nonlinear coefficient of the NLSE are 
derived. But the purpose here is precisely to point out 
explicitly the frequency dependence of the different effec¬ 
tive parameters; and the approximations done in regard 
with it. The propagating equation corresponding to eq. 
1151 in the time domain is: 

d~A(z,t) = iD(id t )A(z,t) +i'y 0 \A(z,t)\ 2 A(z,t) (16) 

Here 70 = 7 (uj re f) is the effective Kerr nonlinearity 
with 7 (w) = n 2 iLo/(cA e ff(u})).(n g /no) 2 . The dispersion 


operator D(id t ) = ^ in> 2 (9”fc)(z9 t ) n /naccounts for dis¬ 
persion at all orders, with t being the retarded time in 
the moving frame at velocity c/n g _ Wo . This NLSE does 
not include any variations of the effective nonlinearity 
over the simulation domain. Effects related to the dis¬ 
persive nonlinearity are then re-introduced using the first 
order perturbation 71 which is accordingly defined as 

71 = dul (w re /). 

d z A(z, t) = iD(id t )A(z, t) 

+ i( r yo + r nid t )\A(z,t)\ 2 A(z,t) (17) 

Eq. [TT] is derived thanks to the fact that the different 
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^/n g (wj = {i j 2 , 3 }) factors in eq. [14] are set outside the in¬ 
tegral sign. Such an operation is only valid if the group 
index n g does not depend on the different frequencies 
uii that contribute to the nonlinear effect. This (math¬ 
ematical) approximation is equivalent to the hypothesis 
according which the strength of nonlinearity only 
depends on the angular frequency which the non¬ 
linear effect takes place at; hence not on the dif¬ 
ferent frequencies contributing to the nonlinear 
effect. There is no evidence that this assumption is valid 
in general. Especially when it comes to nanostructures, 
one may think that the photons that are the most con¬ 
fined contribute the most to the nonlinearity. Although 
the equi-contribution hypo thesis does not hold i n general, 
it can happen that the y/n g (uii)n g (ui 2 )n g (ui 3 ) prefactor 
depends only on uiq. Indeed the different frequencies uii 
are not independent, but must satisfy the energy conser¬ 
vation condition uio + u)\ = W 2 + UJ 3 . In such a case the 
equi-contribution approximation still describes the non¬ 
linear photon dynamics with accuracy; and we will refer 
to such photons as being dispersive photons. 

III. RENORMALIZATION 

As depicted by Fig [T] we shall try to circumvent the 
problem caused by fluctuations of the nonlinear effec¬ 


tive parameters in eq. [14] by finding a proper referential 
wherein the nonlinear response is flat. Looking back at 
eq. [14] tells us that this is indeed the key in order to ob¬ 
tain a time-domain propagation equation like eq. 1161 : the 
different frequencies should have an equi-contribution 
to the nonlinear process; hence no frequency dependent 
pre-factor must appear inside the integral sign of eq. 1141 
To do so, we decided to weight the frequencies by a ad- 
hoc m(ui) contribution. 

Note that we neglect at first any variation of the mode 
field distribution over the pulse bandwidth (gh u o(uii) = 
1 V wo). The reason for this is that variations of the group 
index in PhCWGs account for about 75% of the total 
variation of the effective nonlinear coefficients [H, [ 2 l[ . 
Consequently, dealing with the variations of the slow- 
light factor would be the first step. Moreover this simpli¬ 
fied case is a good test case to present the renormalization 
technique. We will show in next section how the varia¬ 
tions of the mode field distribution can be taken as well 
into account by this technique. 


It appears that a natural choice would be to 
solve the propagation equation for the field ’F(w) = 
yJn g {ui)/noA(ui) so eq. fl4l becomes 


d z ^(z,ui 0 ) = iK(ui 0 )y(z,ui 0 ) +1 ^° n2/ !Va I\H( z ,uh)*^(z, w 2 )4'(z, w 3 )<5 (w2 + w 3 - ui 1 - ui 0 )dui 3 (18) 

cA e f f (ui 0) n 0 J 



FIG. 1. Schematic of the two methods used here to com¬ 
pute the nonlinear pulse propagation in PhC waveguide, a) 
Standard method: the GNLSE is solved for the energy flux, 
b) Input conditions (and the GNLSE) are renormed, and the 
nonlinear propagation is solved for a pseudo-electric 'll field. 

Instead of solving the nonlinear propagation of the 
power flux P = |A| , we are now solving the propagation 
of a pseudo-field *F. Looking in details at the differences 
between eq. [lH] and eq. [15] we see that the slow light en¬ 
hancement factor enhancement S 2 = ( n g /no ) 2 has been 
replaced by S. Interestingly, if we had introduced higher 
order nonlinear effects such as three photons absorption 
(3PA) which are associated to a slow light enhancement 
of S 3 , then the slow light pre-factor would have been 


turned into S as well. A first feature is that the strength 
of the nonlinearity appears to be much weaker. Besides 
if we consider the variations of the Kerr effect 6 ^ 7 ( 0 ;), 
the contribution of the slow-light to the characteristic 
self-steepening time jl 8 } tnl = d ul 'y(ui)/ r y(ui) has been 
halved! Consequently eq. [18] exhibits weaker nonlinear¬ 
ity and even weaker relative nonlinear dispersion than eq. 
M Weaker relative variations of the nonlinearity means 
that the actual nonlinear variations could be taken into 
account with more accuracy by eq. [IS] which could also 
include second order perturbative corrections-. Note that 
although the nonlinearity appears much weaker, the in¬ 
put field has been renormalized as well, hence it is much 
stronger. Consequently global parameters like the soliton 
number are preserved. 

Usually, such a strategy is not convenient because it 
requires additional transformation back and forth be¬ 
tween computed quantities ('F ^ A) and measured ones 
(power flux P = |A| 2 ) -cf. Fig' [U Also one could won¬ 
der whether such a renormalization would not simply 
lead to unphysical solutions, hence simpler equations but 
not describing correctly the physics. Any renormaliza¬ 
tion could be applied to eq. [TT] as long as it is done 
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in consistency with the math -e.g. the integral sign in 
eq. 23; just that not any renormalization would actu¬ 
ally lead to a simpler formulation. Regarding the specific 
choice of m(w) = \fn(w)/no, eq. [9] tells that because 
n g \A\ 2 oc \E \ 2 , then the pseudo-field \F is actually di¬ 
rectly proportional to the Bloch mode electric field. This 
is consistent with the fact that the electric field density 
is indeed the physical quantity that matters for nonlin¬ 
earity, not the power flux. Because the group velocity 
governs most of the physics (slow-light enhancement) in 
PhCWGs, it is then not surprising that it plays a role 
both in the initial normalization of the Bloch field and 
also here in the renormalization of the computed power 
flux. 

In brief, the renormalization technique poses the question 
whether the natural choice that is usually made to com¬ 
pute the evolution of the power flux is right. We think it 
is not for nanostructured systems. 


A. Implementation and comparison 

Nonlinear pulse propagation can now be dealt with in 
two different ways. On one hand, one can use eq. 271 
expanding eventually the Taylor series of the nonlinear¬ 
ity beyond the first order in order to get a better match; 
on the other hand eq. 23 will also provide an accurate 
result, and might be easier to implement. The outcome 
of both equations should be about the same, given that 
they are indeed describing the same physical system. 

To investigate the differences between eq. 23 and eq. 23 
we take as a test case the nonlinear pulse propagation 
experiment performed in conditions similar to those in 
Ref. [16]]. A 2.3 ps Fourier limited pulse is send close to 
the zero group velocity wavelength (ZDW) of an 1.5mm- 
long dispersion-engineered PhC waveguide. The pulse 
peak power is 8 W and corresponds to a soliton number 
of N = 2.1. For consistency with real systems, the equa¬ 
tions have been adapted to include the specificity of op¬ 
tical semiconductor nanostructures, mainly the effect of 
nonlinear absorption and the presence of free carriers fdil - 
[27| . The power flux propagation is computed b y m ean of 
the following generalized Schrodinger eauation|40|. 

^ = -^A-a 3 \A\ 4 A-D(id t )A 
oz 2 

+*(70 + G{idt)) | A\ 2 A — (a + ikoS)N A (19) 

cto = 2 dB/mm and 03 = 25 /W 2 /mm [26] stand 
for the linear propagation loss and three photons ab¬ 
sorption (3PA). The dispersion operator D{idt) = 
En> 2 ffi^)W) n / n accounts for dispersion at all orders, 
with t being the retarded time in the moving frame 
at velocity c/n g (calculated at the input wavelength). 
Besides, we introduce the Kerr operator G{id t ) = 
E n >iffi 7 )W) n / n 'to takes into account the dispersion 
of the Kerr coefficient with the angular pulsation (higher 


order shock terms). Such expansions are intended to pro¬ 
vide accurate numerical results, though it limits the in¬ 
sight into the physical parameters that govern the pulse 
propagation, a and <5 account for the free carriers ab¬ 
sorption and dispersion respectively. Owing to the fact 
that we are considering here a high bandgap material like 
GalnP that exhibits solely three photons (ThPA) and no 
two photons absorption, the self generated plasma does 
not impact much the overall dynamics. Practical details 
related to the way the free carriers effects are computed 
and added to the NLSE are found in Ref. [4p| - [42| . 





Wavelength [nm] Wavelength [nm] 

FIG. 2. Parameters used for the simulation, a) group index 
n g as function of the wavelength, b) Effective Kerr coefficient, 
c) Effective three photons absorption (3PA) coefficient. The 
cyan square indicates the position of the input wavelength in 

Fig E] 

As a guideline we have P 2 = —6.7 ps 2 /mm, /?3 = 
— 1.7 ps 3 /mm, 70 = 2200 /W/m and 71/70 = —170/s, the 
waveguide dispersion and the variation of the nonlinear 
coefficients (Kerr and ThPA) with the angular frequency 
are shown in Fig [2] 

Briefly, the renormed NLSE equation is obtained from 
eq. 23 by (i) expressing it in the frequency domain, (ii) 
dividing the nonlinear coefficient by m(co) n where n = 
2 for effect (e.g. Kerr) and n = 4 for (e.g. 
Three photons absorption), (iii) the input field 'F is then 
obtained by multiplying the initial input field A(uj) by 

Comparison between the results of eq. 23 (red) and its 
counterpart solved for the pseudo electric field \F (blue) 
are shown in Fig[3l The input soliton has undergone a 
blue Soliton-Self Frequency Shift (SSFS) of a few nm and 
dispersive waves are generated in the normal dispersion 
region [ 43 L ItI . In the temporal domain the complex in¬ 
terplay I18i] of the dispersive nonlinearity [45| , free carri¬ 
ers effects [42j, and the SSFS, results in a pulse advance 
of a few picoseconds. We see the spectral position of 
the dispersive wave (in the normal dispersion region) and 
the time advance (lOps vs. 15ps) differ between the two 
models. Indeed after the initial generation, the dispersive 
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AX [nm] 



FIG. 3. Comparison between the two numerical models, a) 
Spectra, b) Temporal trace, dashed-black: input pulse; blue: 
with renormalization m = \Jn g {uo)/no- red: without renor¬ 
malization. Ao indicates the position of the zero group veloc¬ 
ity dispersion point (ZVD). 


wave interacts with the soliton through cross-phase mod¬ 
ulation and cascaded four wave mixing (461]: therefore 
its amplitude and position depends greatly on the exact 
form of the Kerr nonlinearity. The amplitude of the dis¬ 
persive waves grows much stronger when the renormal¬ 
ization procedure is employed; and because the spectral 
recoil appears in reaction to the emission of the disper¬ 
sive wave, the SSFS is much stronger for the renormalized 
case. Divergences between the two numerical models are 
clearly visible. 


B. Comparison with an analytic solution 

Now that we have shown that the two approaches (the 
nominal GNLSE and its renormalized counterpart) lead 
to different results, we must determine which model is 
(the most) correct. The case that we just discussed cor¬ 
responds to a realistic case where both the medium (dis¬ 
persion, nonlinearity, absorption) and the input parame¬ 
ters (power, duration) are within reach of current exper¬ 
iments. So it would be in principle possible to perform 
such an experiment and compare the two models with 
the experimental results. However such measurements 
are not available yet. Another way would be to confront 
directly the results of the two models with an analytic 
case. 

Recent work by Erkintalo et al. [46| demonstrated that 
nonlinear pulse propagation (continuous spectrum) and 
cascaded FWM (discrete spectrum) are closely related 
and that the nonlinear pulse dynamics can be described 
as a cascade of FWM events. This means that the ca¬ 
pacity of an equation to accurately reproduce the reality 
is intrinsic to its capacity to deal correctly with FWM. 
Although no analytic solution exists for the exact case we 
just studied we can still simplify our problem to the pro¬ 


pagation in a lossless and L = 100gm short PhCWG of 
a ?o = 8W single continuous wave beam of same central 
frequency ( oo re f ) as previously. In such a situation the 
short length of the waveguide and the absence of prop¬ 
agation loss render the undepleted pump approximation 
valid. The FWM conversion efficiency [7j depending on 
the pump signal detuning (<5w) is then expressed through 
the analytic formulation: 


V(6w) = ( 7 FWM{5u)P 0 L)\ Sinh \f“ )L) ) 2 (20) 

g(ou})L 

g 2 (5uj) = (7 fwm{Suj)P 0 ) 2 

— (AK L uRef(Suj)+AK NLuiRe f(5uj)) 2 /4: 

( 21 ) 

AKLuRef (^^) = 2AT(cu re y) A (cu re j T Jcu) 

— K{uj re f — Slo) 

( 22 ) 


AK NLu ,R e f(5u 1 ) = 2P 0 ('Jxpm(Suj) + 1xpm(—Suj) 

- 1spm{u re /)) (23) 


7 fwm, lx pm, is pm account respectively for the 
FWM nonlinear coupling coefficient, the cross phase 
modulation (XPM) between the strong pump and the 
weak signal/idler and the self-phase modulation (SPM) 
of the pump [2]J. Especially it takes into account the 
different modes fields overlap. According to the notation 
used previously, these effective nonlinear coefficients are 
expressed as: 


IS pm (w re /)= r. e/ : 2J .^ cure fiuiref, ^ref, ^ref ) 

C -Aef f\Uref) 


( n g(u)ref ) x2 


7 xpm{w)= 


n 0 
tori2i 


(24) 


cA e ff(uj) 
n g (u> re f)n g (uj) 


gh 

tore f(uiref, LJ re f) 

) (25) 


/ . , (Wre f+Suj)n 2 i u , £. ^ 

T FWM T 7 7 d^LOre t \^ref OCU, UJ re f , Cdref) 

cA e ff(CU+dLU) J 


Tig (id re f ) \J Tlg{u) re f 8bj)flg {id r e / -\~duj ) 

V T~2 ) 


(26) 


These coefficients are computed for each pump-signal 
detuning 5w\ as a result we get the FWM gain curve 
as shown in thick black in Fig [I] We compare now this 
analytic curve to the results given by the different models. 


First we see that the GNLSE (thick dashed green) does 
not appear to converge any better than the standard 
NLSE (light dashed blue). The NLSE does not take 
into account any variations of the nonlinearity, hence 
IS pm = lx pm = 1 fwm- While the NLSE tends to 















FIG. 4. r? = Pidier^L)/P aigrl ai{ 0) for various pump-signal de¬ 
tuning after a propagation of L = 100/rm. Black: Analytic 
model. Thin dashed blue: results of the NLSE. Thick dashed 
green: GNLSE. red: GNLSE renormalized where only the 
group index variations are taken into account, cyan dots: 
GNLSE where both slow light and effective area variations 
are included into the renormalization. 


overestimate the FWM bandwidth, the GNLSE under¬ 
estimates the FWM gain. Thus the inclusion of self- 
steepening (i.e. the dispersion of the nonlinearity) does 
not appear as a great improvement. One must note that 
the GNLSE is still essential to model effect specifically 
related to self-steepening like the formation of a shock 
front or to explain the e nerg y dependent time-advance of 
the nonlinear pulse (Tl, |45| . We see that the renormal¬ 
ization of the slow light enhancement factor (thick red) 
improves the convergence of the GNLSE: for small de¬ 
tuning (<5A < 2 nm) the GNLSE and the analytic model 
converge; for larger detuning some discrepancies still re¬ 
main, but the overall error is still weaker than for the 
un-renormalized GNLSE or the simple NLSE. The re¬ 
maining error is due to the fact that the effect of the 
variation of the effective modal area is neglected at first. 
We see that if we then renormalized the GNLSE to take 
as well into account that later effect (cyan dots), we ob¬ 
tain a very good agreement between the analytic model 
and the renormalized GNLSE. We present in the last 
section how the renormalization technique is general and 
can take into account variation of the modal area. 

This demonstrates that the renormalization technique 
that we present here really constitutes an improvement 
compared to previous formulations. More generally our 
technique could also be seen as a generalization of Ref. 


21| | that deals in a seamless way with the generation of 


multiple signal and idler orders, hence not only limited 
to a discrete set of a few (usually four) beams. 


C. Discussion 

The differences between the two models (with or with¬ 
out the renormalization) indicate that the physics gov¬ 
erning the nonlinear pulse is necessarily different. Espe¬ 
cially, one of the largest changes between eq. [T7] and eq. 


[T8l lies in the ratio of Self-Phase (SPM) to Cross-Phase 
Modulation (XPM) intensity. Usually the XPM is twice 
as many times stronger as the SPM for a given frequency. 
After the renormalization of the slow light variations, the 
photons are weighted by the quantity m(oj) 2 = n g (uj)/no 
and they do not have the same contribution to the nonli¬ 
near index change. The weight factor corresponds to the 
increase of the photonic density of states in the PhCWGs 
compared to the bulk material. Indeed the first deriva¬ 
tive of the band diagram could be both interpreted as the 
group velocity or the density of states: photons with a 
higher density of states will have a higher probability to 
interact with the nonlinear medium. By analogy, the in¬ 
clusion of the variations of the effective modal area would 
correspond to dividing the photon eigen contribution by 
their effective volume. Thus the enhancement of the non¬ 
linearity due to the slow light and the tight confinement 
of light can be seen as an enhancement caused by the in¬ 
crease of the optical density of states. The photons with a 
high density of states are also less influenced by the other 
photons with a lower density of states ( SPM < 2 XPM). 

Another interesting effect lies in the magnitude of the 
shock term : the stronger this term is, the faster will 
the pulse form a shock front. The formation of an opti¬ 
cal shock -more precisely the presence of the Kerr shock 
term- might play a very important role in the genera¬ 
tion of dispersive linear waves (DSW) that could then be 
generated without strong dispersion requirements, for in¬ 
stance the presence of ZVD point is not mandatory [47j. 
The intensity of the XPM and its dependence with the 
angular frequency plays then a predominant role. After 
the renormalization, the dispersive nonlinearity ( 71 / 70 ) 
in PhCWGs is only about one third of its nominal value 
in eq. HU while at the same time the weight of the pho¬ 
tons has been strongly modified. Consequently the be¬ 
havior of PhCWGs with regards to this new DSW gen¬ 
eration scenario is different to what is found in other 
systems. More generally, the renormalization redefines 
in a non-trivial way the interaction, mediated by the ma¬ 
terial nonlinearity, between the photons. Consequently, 
in the laboratory reference frame (i.e. considering only 
the power flux and no renormalization) the PhCWG be¬ 
havior would be different to what is a priori expected. 


IV. IMPACT OF MODE AREA 

Thus far, we only presented how to include the varia¬ 
tions of the slow light factor in the renormalization pro¬ 
cess. The main reason is that the slow-light is responsible 
for most of nonlinearity variations in PhCWGs. However 
the change in the effective modal area still accounts for 
about 25% of the dispersive nonlinearity; and we have 
shown through the comparison with an analytic test case 
that it has a noticeable impact on the overall nonlinear 
dynamics. Consequently we now show how the renor¬ 
malization method is also able to include variations of 
the effective modal area. 
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By analogy with what has been done in the previous 
section, the renormalization procedure is applicable as 
well for the modal area subject that gh u o can be decom¬ 
posed as ghuo(wi,W 2 ,W 3 ) = g{ui0)h(oji)h(oj2)h(uj 3 ). If 


such is the case, then we define the renormalization func¬ 
tion m(w) = h{uj)^n g ( u> ), and the propagation equation 
to solve becomes: 


d z 'S>(z,uio) = iK(uj 0 )^(z,ujo) +i 


LUon 2 i n 


guio 


cA e ff(ui 0) no 


g(ujo)h(LUo) / ^(z, 003)6(1^2 + W3 - wi - uo)duf 

(27) 


Consequently h(u>) weights the individual photons’ 
contribution to the Kerr nonlinearity: it includes implic¬ 
itly most of the dispersive nonlinearity (in PhCWGs). 
On the contrary g(cuo)h(u>o)/A e ff (u>o ) stands for the 
dispersive part of the modal area. The problem of 
the factor decomposition of gh u 0 is directly related to 
the question whether the nonlinearity in PhCWGs can 
be modelled accurately using an analytic function j20j . 
It has been demonstrated that, for dispersion engineered 
PhCWGs like the one considered in the present paper, 
the nonlinearity can be fitted by a Morse type potential 
function with four adjustable parameters. However the 
decomposition of such function in factor decomposition 
will only gives an approximate. More generally the best 
way to decompose the variations of the effective modal 
area is still an open question. 

In any case, it is always possible to choose a de¬ 
composition that preserves the self-phase modulation 
(. SPM{uj ) oc g(ui)h(oj) 3 ) and cross-phase modulation 
XPM(co re f,ui) oc g(uj)h(iO re f)' 2 h(io). oj re f is defined as 
the center of the frequency domain. 
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FIG. 5. a) Color map: h(co) computed for different u> re f (y- 
axis). b) h(uj) averaged on the different a ] re f for the normal 
dispersion region (blue); and the anomalous dispersion region 
(red). Shaded area indicates the standard deviation of the 
average, c) Same as b), but showing the dispersive part of 
the nonlinearity g{uj)h(ui) / A e ff(uo). 


h{u) 

/ SPM(ui) 

(28) 

I XPM (u re f, LO) 

9 M =4 

XPM(tUref,Uj) 3 

1 SPM(ui) 3 

(29) 


If the decomposition is consistent -i.e. if ghu>i can in¬ 
deed be decomposed in factors-, h(u) and g(u) do not 
depend on the central frequency u re f that is chosen. 
Otherwise, only the XPM created by a pulse centered at 
ui re f, as well as the SPM for any frequency, are included 
correctly. Such approximation could still be sufficient if 
the propagation is dominated by a single strong pulse. 
Checking how h{ui re f) and g(u} re f) depends on w re / is 
essential to assess the validity of the renormalization for 
taking into account effects related to the effective modal 
area. In Fig[5}a), we show the value of h(ui) computed 
according to eq. [2H1 depending on the central frequency 
uiref (Y-axis). 

We observe two main zones: one ranging 1525 — 
1560nm (anomalous dispersion) and the other one rang¬ 
ing 1560 — 1610nm (ranging from the first zero GVD to 


the band edge). The two zones are separated by the first 
zero group velocity dispersion (GVD) point (cf. Fig [2} 
a)). Inside each zone, g(ui)h(co) and h(u>) do not depend 
on the central reference frequency as seen in Fig[5]-b,c). 
This brings forth two major conclusions. First it ap¬ 
pears that inside each zone, the decomposition of gh(u>i) 
holds; and it is therefore possible to describe correctly 
the pulse evolution using the renormalization technique 
(of course providing that the simulation spectral domain 
being confined within one zone). Second, the presence of 
two distinct and well defined zones indicates that there 
is actually a change in the physics governing the photons 
evolution. 

For high frequencies (small wavelengths), the photons 
have a nonlinear dispersive behavior corresponding to the 
fact that h(ui)g(ui)/A e ff(uj ) varies while h(u>) is almost 
constant (Fig[5]-b), blue). This indicates that the photons 
have a quasi equi-contribution to the nonlinearity. On 
the contrary, for low frequencies (long wavelengths), the 
individual contribution of photons is more pronounced 
(FigEJb), red) and h(u>) exhibits variations up to 60% 
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while h(uj)g(uj) remains flat. Weighting the eigen pho¬ 
tons contribution to the nonlinearity is essential. 

In dispersion-engineered PhCWGs, the slow light at long 
wavelength (low frequency) is caused by the presence 
of a complete Photonic Bandgap (PBG). Close to the 
band edge, the physics is then similar to what is found in 
Bragg Gratings j48|. In contrast, the slow-light obtained 
at higher frequency through dispersion-engineering has 
another nature and arises thanks to the complex inter¬ 
ferences occurring inside the Bloch mode. Although the 
decomposition of gh(cu) into a product of functions is not 
mathematically exact, we see here that it could never¬ 
theless be a useful metric to sort-out slow light [49[ into 
categories depending on the equi-contribution to the non¬ 
linearity that the photons have -or do not have-. 

In the introduction, we disregarded the coupled set 
of equations as a possible solution for modeling systems 
with large variation of the nonlinear parameters. The 
fact that the photons are split here into two well defined 
domains tends to rehabilitate a posteriori such strategy. 
Still one must be careful on the way photons are taken 
into account at the limit between the two domains, es¬ 
pecially considering that this point is precisely the zero 
group velocity dispersion (ZVD) point. 

Finally we focused on SPM, XPM in the factor decom¬ 
position of gh(u>i). The photon weight h(uj) depends on 
this peculiar choice. Therefore a material exhibiting a 
different nonlinear response, like the presence of a 
nonlinearity or a different form for the tensor, would 
lead to a different /i(w). 


V. CONCLUSION 

We have presented a new method to incorporate in an 
efficient way into the GNLSE the variations of the nonlin¬ 
earity that exist in systems with structured slow-light: 
effects linked to variation of the slow-light enhancement 
factor could be taken into account through renormaliza¬ 
tion. Fluctuations of the mode effective area can be dealt 
with by this technique as well. This would be of impor¬ 
tance especially for system like nanowires, where vari¬ 
ations of the effective mode area are important. As a 
result this paper gives practical hints regarding the way 
nonlinear pulse propagation in nanostructured systems 
could be computed; and what could be the limitation of 
current models. 

This is crucial as more studies are precisely focusing on 


higher order nonlinear effects and their mutual interplay 
flq . flit l42l . l45l H3 . Our model is consistent with an an¬ 
alytic set of equations derived to model discrete FWM 
events [2lj|; and could be considered as a generalization 
of that article. Besides, it is worth to note that our tech¬ 
nique does not increase the computational burden com¬ 
pared to the resolution of a standard GNLSE. 

This study was also the occasion for more fundamen¬ 
tal considerations. In particular we found relevant not to 
compute directly the propagation of the power flux but 
to weight first the photon contribution by their optical 
density of states. Although the modal area depends on 
the considered nonlinear effect and its associated tensor, 
we found two classes of slow light in dispersion-engineered 
PhCWGs: engineered slow light with a normal nonlinear 
dispersive behavior, and a region close to the photonic 
band gap where the weight factor of the photons con¬ 
tributes greatly. Such study might change how we per¬ 
ceive and understand slow-light effects which are in fact 
more related to high density of states light effect. 

Usually, the improvement of the models which deal 
with nonlinear pulse propagation comes along with the 
addition of extra operators in order to describe the new 
effects. Here the renormalized equation keeps the for¬ 
mulation with no extra terms added. Indeed, the new 
phenomena that are observed lie inside the renormaliza¬ 
tion function m(oj ), not in the GNLSE itself. Namely 
the unbalance between SPM and XPM could be in¬ 
terpreted as inertial forces that appear because of the 
non-trivial relationship between the laboratory referen¬ 
tial and the PhCWGs one, where some photons appear 
more immune to perturbation or prone to perturb others. 
Within this new reference frame, all the semi-analytic 
method/models developed so far, like the momentum 
method, remain valid. 

Finally the present discussion only focused on 
PhCWGs, a system where the nonlinear variations are 
extreme. Part of our conclusions would anyway also 
apply to other nano-structured systems like nanowires 
where the relatively weaker nonlinear variations must 
be seen in regards to the very large optical bandwidth 
these systems support. Especially we have shown that 
variations of the effective area -which are dominant over 
slow-light in nanowires- can be taken into account by 
the renormalization method; and that this method can 
describe with accuracy both self-phase and cross phase 
modulation effects. 
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